#include <string>
#include <algorithm>
#include <numeric>
#include <istream>
#include <vector>
#include <fstream>
#include <deque>
#include <cassert>
#include <map>
#include <set>
#include <stdexcept>
#include <ext/hash_map>
#include <ext/hash_set>


// boost functions
#include <boost/lexical_cast.hpp>
#include <boost/filesystem/operations.hpp>

// my functions
#include <portable_timer.h>
#include <progress_indicator.h>
#include <blast_results.h>
#include <stlhelper.h>
#include <print_error.h>
#include <ultoa.h>
#include <index_string.h>
#include <get_filename.h>
#include <single_linkage.h>
#include <open_fstream.h>


#include "get_arg.h"

#include <gc.h>
using std::make_pair;
using std::map;
using std::cerr;
using std::ifstream;
using std::ofstream;
using std::pair;
using std::cout;
using std::string;
using std::vector;
using std::deque;
using std::sort;
using boost::lexical_cast;
using std::ostringstream;
using stlhelper::lookup;




struct  t_prot_pair_blast_res
{
	unsigned prot_id1;
	unsigned prot_id2;
	float	 bitscore;
	float	 e__value;
	unsigned align_len;
	t_prot_pair_blast_res(	unsigned prot_id1_,
							unsigned prot_id2_,
							t_blast_results results):
							prot_id1( prot_id1_),
							prot_id2( prot_id2_),
							bitscore( results.bitscore),
							e__value( results.e__value),
							align_len(results.align_len){}
};

bool operator < (const t_prot_pair_blast_res& a, const t_prot_pair_blast_res& b)
{
	if (a.bitscore < b.bitscore)
		return true;

	// if same evalue
	else if (a.bitscore == a.bitscore)
	{
		//only save if lower e_value
		if (a.e__value > b.e__value)
			return true;
		else if (a.e__value == b.e__value)
		{
			//  or shorter length
			if (a.align_len > b.align_len)
				return true;
		}
	}
	return false;
}

typedef __gnu_cxx::hash_map<unsigned, unsigned>	t_blast_index_per_gene;
typedef __gnu_cxx::hash_set<unsigned>			t_index_set;
typedef deque<t_blast_index_per_gene >			t_matrix_blast_indices;
typedef t_blast_index_per_gene::iterator		t_blast_index_per_gene_iterator;
typedef t_blast_index_per_gene::const_iterator	t_blast_index_per_gene_citerator;
typedef deque<t_prot_pair_blast_res>			t_blast_res_array;
//typedef deque<std::map<unsigned, unsigned> >		t_matrix_blast_indices;
//std::map<unsigned, unsigned>::iterator		t_blast_res_index_iterator;



//________________________________________________________________________________________

//		parse_protid_gene_id_from_file

//________________________________________________________________________________________
void parse_protid_gene_id_from_file(const string&			file_name,
									deque<unsigned>&		prot_id_to_gene_id,
									t_index_string&			prot_ids,
									t_index_string&			gene_ids,
									deque<char>&			gene_ids_single_prot,
									t_program_args&			args)
{
	ifstream is;
	open_fstream(is, file_name);

	unsigned cnt_lines = 0;

	// get data for one gene
	string prot_name, gene_name;
	is >> prot_name >> gene_name;
	is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
	args.VERBOSE() << "\tReading gene to protein mapping...\n";
	t_progress_indicator 				dots(args.VERBOSE(),  20000);
	while (is)
	{
		dots++;
		cnt_lines++;
		bool added;
		unsigned prot_id = prot_ids.index(prot_name, added);
		if (added)
		{
			unsigned gene_id = gene_ids.index(gene_name, added);

			// has a transcript for this gene been encountered before.
			// if not, it may be a single transcript gene
			if (added)
			{
				assert(gene_id == gene_ids_single_prot.size());
				gene_ids_single_prot.push_back(1);
			}
			else
				gene_ids_single_prot[gene_id] = 0;

			// save gene_id lookup
			prot_id_to_gene_id.push_back(gene_id);

		}
		else
		{
			args.CERR() << "WARNING:\tDuplicated prot ID " << prot_name
						<< " with previous gene of " << gene_ids.string_at(prot_id_to_gene_id[prot_id])
						<< " and now with gene of " << gene_name << "\n";
		}

		// get next line
		is >> prot_name >> gene_name;
		is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
	}

	if (is.eof())
	{
		dots.finish();
		args.VERBOSE(2) << "\tFinished parsing all possible identifiers:\n";
		args.VERBOSE(2) << "\t\t" << prot_ids.size() << " possible protein transcript identifiers.\n";
		args.VERBOSE(2) << "\t\t" << gene_ids.size() << " possible gene identifiers, of these \n";
		assert(gene_ids.size() == gene_ids_single_prot.size());
		unsigned cnt_single_transcript_gene =
			std::accumulate(gene_ids_single_prot.begin(),  gene_ids_single_prot.end(), 0U);
		args.VERBOSE(2) << "\t\t" << cnt_single_transcript_gene << " are for genes with single transcripts.\n";
		return;
	}

	throw std::runtime_error("Invalid protein / gene name on line #" +
								ultoa(cnt_lines) + " of " + file_name);


}
















//________________________________________________________________________________________

//		lookup_prot_name

//________________________________________________________________________________________
inline unsigned lookup_prot_name(t_index_string& prot_ids, const string& prot_name)
{
	unsigned prot_id;
	if (!prot_ids.is_indexed(prot_name, prot_id))
		throw std::runtime_error("Protein identifier in blast result [" + prot_name +
								"] not found in the supplied protein id to gene id "
								"lookup table.");
	return prot_id;
}

//________________________________________________________________________________________

//		save_best_hitpair_for_multiple_transcripts

//________________________________________________________________________________________
void save_best_hitpair_for_multiple_transcripts(t_blast_res_array&  			blast_res,
												t_matrix_blast_indices&			matrix_blast_indices,
												unsigned	  					gene_id1,
												unsigned	  		 			gene_id2,
												const t_prot_pair_blast_res&	new_blast_res)
{
	// new blast result for gene_id1 vs gene_id2
	//		save it in blast_res
	//		save index of blast_res in matrix_blast_indices[gene_id1][gene_id2]
	if (matrix_blast_indices[gene_id1].count(gene_id2) == 0)
	{
		// save index into new blast results in index (save index both ways so can look up
		matrix_blast_indices[gene_id1][gene_id2] =
		matrix_blast_indices[gene_id2][gene_id1] = blast_res.size();
		blast_res.push_back(new_blast_res);
		assert(blast_res[matrix_blast_indices[gene_id1][gene_id2]].prot_id1 ==
					new_blast_res.prot_id1);
		assert(blast_res[matrix_blast_indices[gene_id1][gene_id2]].prot_id2 ==
					new_blast_res.prot_id2);
		return;
	}

	// There is already a blast result for gene_id1 vs gene_id2
	//		retrieve it from blast_res (index == matrix_blast_indices[gene_id1][gene_id2])
	// Save it if better
	unsigned blast_res_index = matrix_blast_indices[gene_id1][gene_id2];
	assert(blast_res_index < blast_res.size());
	if (blast_res[blast_res_index] < new_blast_res)
		blast_res[blast_res_index] = new_blast_res;
}


//________________________________________________________________________________________

//		print_saved_pairs_with_multi_transcripts

//________________________________________________________________________________________
void print_saved_pairs_with_multi_transcripts(	const t_index_string&			prot_ids,
												const t_index_string&			gene_ids,
												const deque<unsigned>&			prot_id_to_gene_id,
												const t_blast_res_array&		blast_res,
												const t_matrix_blast_indices&	matrix_blast_indices,
												deque<char>&					gene_ids_unprinted,
												t_program_args&					args)
{
	args.VERBOSE() << "\tPrinting best hitting transcripts for multiply spliced genes...\n";
	unsigned DOTS_SIZE = 500000;
	t_progress_indicator dots(args.VERBOSE(),  DOTS_SIZE);
	dots.use_timer();
	for (unsigned gene_id1 = 0; gene_id1 < matrix_blast_indices.size(); ++gene_id1)
	{
		t_matrix_blast_indices::value_type::const_iterator jj = matrix_blast_indices[gene_id1].begin();
		t_matrix_blast_indices::value_type::const_iterator end = matrix_blast_indices[gene_id1].end();
		for (; jj != end; ++jj)
		{
			unsigned gene_id2 = jj->first;
			unsigned blast_index = jj->second;
			if (gene_id2 > gene_id1)
				continue;
			++dots;
			const t_prot_pair_blast_res& res = blast_res[blast_index];
			cout	<< prot_ids.string_at(res.prot_id1) << "\t"
					<< prot_ids.string_at(res.prot_id2) << "\n";
			gene_ids_unprinted[gene_id1] = 
			gene_ids_unprinted[gene_id2] = 0;
		}
	}

	// count and return all genes unprinted
	{
		unsigned cnt_unprinted_genes_left =
			std::accumulate(gene_ids_unprinted.begin(),  gene_ids_unprinted.end(), 0U);
		if (cnt_unprinted_genes_left)
			args.VERBOSE(2) << "\t\t" << cnt_unprinted_genes_left << " genes have not been printed.\n";
		else
			args.VERBOSE(2) << "\t\tAll genes have been printed!\n";
	}	
}




//________________________________________________________________________________________

//		print_matching

//________________________________________________________________________________________
void parse_blast_res_from_file(	const string&			file_name,
								deque<unsigned>&		prot_id_to_gene_id,
								t_index_string&			prot_ids,
								t_index_string&			gene_ids,

								// multiple transcript genes saved
								t_blast_res_array&		multi_transcripts_blast_res,
								t_matrix_blast_indices& matrix_multi_transcripts_blast_indices,

								// single transcript genes printed out immediately
								deque<char>&			gene_ids_single_prot,
								t_program_args&			args)
{
	args.VERBOSE() << "\tParsing blast results from [" << 
					get_filename(file_name) << "]...\n";

	ifstream is;
	unsigned file_size = open_fstream(is, file_name);

	unsigned DOTS_SIZE = 50000;
	t_progress_indicator dots(args.VERBOSE(),  DOTS_SIZE);
	dots.use_timer();
	
	unsigned	cnt_lines = 0;
	unsigned	cnt_single_transcript_gene_results = 0;
	unsigned	cnt_multi_transcript_gene_results = 0;

	deque<char>		gene_ids_single_prot_printed(gene_ids.size(), 0);
	deque<char>		gene_ids_in_blast_res(gene_ids.size(), 0);
	try
	{
		string line;
		while (is)
		{
			++cnt_lines;

			//
			// update dots
			//
			++dots;
			if (dots.value() % (DOTS_SIZE * 10) == 0)
			{
				boost::intmax_t file_size_so_far = is.tellg();
				unsigned dots_needed = unsigned(dots.value() * 1.0 *
												file_size / file_size_so_far);
				dots_needed += DOTS_SIZE;
				dots.set_limit(dots_needed);
			}
	
			// get data for one protein pair
			string prot_name1, prot_name2;
			is >> prot_name1 >> prot_name2;
			t_blast_results new_blast_res;
			is >> new_blast_res;
	
			if (!is.good())
				break;
	
	
			// get protein and gene indices
			unsigned prot_id1 = lookup_prot_name(prot_ids, prot_name1);
			unsigned prot_id2 = lookup_prot_name(prot_ids, prot_name2);
			unsigned gene_id1 = prot_id_to_gene_id[prot_id1];
			unsigned gene_id2 = prot_id_to_gene_id[prot_id2];
	
	
			// ignore comparisons between transcripts of the same gene
			if (gene_id1 == gene_id2)
				continue;

			gene_ids_in_blast_res[gene_id1] = 1;
			gene_ids_in_blast_res[gene_id2] = 1;

			// sort gene indices so that comparisons don't have to be symmetrical
			if (gene_id1 > gene_id2)
			{
				std::swap(gene_id1, gene_id2);
				std::swap(prot_id1, prot_id2);
			}

			// no more transcripts for this gene, just print out
			if (gene_ids_single_prot[gene_id1] && gene_ids_single_prot[gene_id2])
			{
				cout
					<< prot_name1 << "\t"
					<< prot_name2 << "\n";
				++cnt_single_transcript_gene_results;
				gene_ids_single_prot_printed[gene_id1] = 1;
				gene_ids_single_prot_printed[gene_id2] = 1;
			}
			else
			// save if necessary
			{
				save_best_hitpair_for_multiple_transcripts(	multi_transcripts_blast_res,
															matrix_multi_transcripts_blast_indices,
															gene_id1,
															gene_id2,
															t_prot_pair_blast_res(	prot_id1,
																					prot_id2,
																					new_blast_res));
				++cnt_multi_transcript_gene_results;
			}

			
		}
		dots.finish();
	}
	catch(std::runtime_error& e)
	{
		throw std::runtime_error("Invalid blast results format on line #" +
									ultoa(cnt_lines) + " of " + file_name + "\n" + e.what());
	}
	if (!is.eof())
		throw std::runtime_error("Could not read from the file " + file_name);
	args.VERBOSE(2) << "\t" << string(60, '_') << "\n";
	args.VERBOSE(2) << "\t" << cnt_lines << " blast results read from [" << file_name << "].\n";
	args.VERBOSE() << "\t" << cnt_single_transcript_gene_results	<< " blast results from genes with single transcripts.\n";
	args.VERBOSE() << "\t" << cnt_multi_transcript_gene_results		<< " blast results for genes with multiple transcript.\n";

	// count gene_ids in blast
	{
		unsigned cnt_gene_ids_in_blast_res = 
			std::accumulate(gene_ids_in_blast_res.begin(),  gene_ids_in_blast_res.end(), 0U);
		args.VERBOSE(2) << "\t\t" << cnt_gene_ids_in_blast_res << " genes in blast result pairs.\n";
	}
	// count single transcript printed
	{
		unsigned cnt_gene_ids_in_blast_res = 
			std::accumulate(gene_ids_single_prot_printed.begin(),  gene_ids_single_prot_printed.end(), 0U);
		args.VERBOSE(2) << "\t\t" << cnt_gene_ids_in_blast_res << " single transcript genes printed.\n";
	}

	// count and return all genes unprinted
	{
		for (unsigned i = 0; i < gene_ids_single_prot.size(); ++i)
			gene_ids_single_prot[i] = gene_ids_in_blast_res[i] - gene_ids_single_prot_printed[i];
		unsigned cnt_unprinted_genes_left =
			std::accumulate(gene_ids_single_prot.begin(),  gene_ids_single_prot.end(), 0U);
		args.VERBOSE(2) << "\t\t" << cnt_unprinted_genes_left << " genes remaining to be printed.\n";
	}	

}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_portable_timer timer;

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// get command line arguments
	t_program_args args;
	try
	{
		if (!process_command_line_options(argc, argv, args))
			return 1;
	
		//
		//	protein to gene id lookup
		//
		deque<unsigned>					prot_id_to_gene_id;
		t_index_string					prot_ids;
		t_index_string					gene_ids;
		deque<char>						gene_ids_single_prot;
		parse_protid_gene_id_from_file(	args.f_prot_to_gene_id,
										prot_id_to_gene_id,
										prot_ids,
										gene_ids,
										gene_ids_single_prot,
										args);
	
	
	
		//
		//	parse blast results
		//
		//
		t_blast_res_array		multi_transcripts_blast_res;
		t_matrix_blast_indices	matrix_multi_transcripts_blast_indices(gene_ids.size());
		parse_blast_res_from_file(	args.f_input,
									prot_id_to_gene_id,
									prot_ids,
									gene_ids,
									multi_transcripts_blast_res,
									matrix_multi_transcripts_blast_indices,
									gene_ids_single_prot,
									args);
	
	
		//
		//	parse best hitting transcript from multiply spliced genes
		//
		//
		deque<char>	gene_ids_unprinted;
		gene_ids_unprinted.swap(gene_ids_single_prot);	// change name to make intent clearer
		print_saved_pairs_with_multi_transcripts(prot_ids,
												 gene_ids,
												 prot_id_to_gene_id,
												 multi_transcripts_blast_res,
												 matrix_multi_transcripts_blast_indices,
												 gene_ids_unprinted,
												 args);
	
	
	
		args.VERBOSE() << "\tAll finished after " << timer.pretty_reset() << "\n";
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}



